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multi-dimensional feedback delay 
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Abstract 


In this paper, we present a new method for solving fractional optimal 
control problems with delays in state and control. This method is based 
upon Bernstein polynomials basis and feedback control. The main advantage 
of feedback or closed-loop control is that one can monitor the effect of such 
control on the system and modify the output accordingly. In this work, we 
use Bernstein polynomials to transform the fractional time-varying multi- 
dimensional optimal control system with both state and control delays, into 
an algabric system in terms of the Bernstein coefficients approximating state 
and control functions. We use Caputo derivative of degree 0 < a < 1 as the 
fractional derivative in our work. Finally, some numerical examples are given 
to illustrate the effectiveness of this method. 


Keywords: Delay fractional optimal control problem; Caputo fractional 
derivative; Bernstein polynomial. 


1 Introduction 


The general definition of an optimal control problem requires the minimiza- 
tion of a functional over an admissible set of control and state functions sub- 
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ject to dynamic constraints on the states and controls. A Fractional Optimal 
Control Problem (FOCP) is an optimal control problem in which either the 
performance index or the differential equations governing the dynamic of the 
system or both contain at least one fractional order derivative term [1, 2, 17]. 


Fractional Differential Equations ( FDEs ) have been the focus of many 
studies due to their appearance in various applications in real-world physical 
systems. For example, it has been illustrated that materials with memory 
and hereditary effects and dynamical processes including gas diffusion and 
heat conduction can be more adequately modeled by FDEs than integer- 
order differential equations [13, 18, 20]. Some other applications of FDEs 
are in behaviors of viscoelastic materials, biomechanics and electrochemical 
processes ( see [3, 5] for more details ). 


Most FOCPs do not have exact solutions, so in these cases approximation 
methods and numerical techniques must be used. Recently, several approxi- 
mation methods to solve FOCPs have been introduced [4, 14, 18]. 


Real life phenomena have been described more precisely by Delay Differ- 
ential Equations, so Delay Fractional Optimal Control Problem ( DFOCP ) 
has become the focus of many researchers in the last decade. Baleanu in [6] 
and Jarad in [11] analyzed the fractional variational principles for some kinds 
of DFOCPs within Riemann-Liouville and Caputo fractional derivatives re- 
spectively and made their corresponding Euler-Lagrange equations. In this 
paper, we present a novel strategy based on Bernstein polynomials (BPs) to 
solve DFOCPs. Consider the following DFOCP 


Min J=3 fo [2TH Q(é)2(t) + uT (t)R(t)u(t)ldt, (1) 
St 
pDexi(t) = Uh 44,7 (t) xj (t) + Up 1i,n(t)ue(t) (2) 
+25 _ 1 (@a)i,5 (taj (t — m1) + Up (ba)i,n(Qua(t— 2), l<Sisr, 
£;(t) = 250, teE [—m1, 0], 1<j<r, 3 
ux(t) = uro, t € [—n2, 0}, 1l<k<s, ( ) 


where x(t) = [x1(t)---a,(t)|? and u(t) = [ui(t)---us(t)]” are respectively 
the state and control functions. Also, Q(t) and R(t) are respectively, r x r 
and s X s semi-positive and positive definite time-varying matrices of the 
state and control’s coefficients in the cost function with continuous functions 
as their entries. Furthermore, a;,;(t), (@a)i,j;(t), bi,n(t) and (ba)i,4(t) are 
continuous functions which are respectively the coefficients of 7; (t), x;(t — 
m) for (1 < 7 < r) and uxg(t), ux(t — n2) for (1 < k < 8s) in the ith 
fractional differential equation (2) and 7,72 > 0 are given constant delays. 
The fractional derivative is defined in Caputo sense, i.e. 


dr 


t —a 
Mia) Je (t 1) “4a;(r)dr, 0<a <1, 
Li, a=. 


oD ai(t) = (4) 
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In the numerical solution of dynamical systems, polynomials or piecewise 
polynomial functions are often used to present the approximate solutions [9, 
10, 21]. The effectiveness of using Bernstein polynomials for solving FOCPs 
have been demonstrated before [4, 14]. In the present paper, we seek an 
optimal feedback control function to find the approximate solution of DFOCP 
(1) - (3) by using Bernstein polynomials. 

This paper is organized as follows. In Section 2 we give some preliminiaries 
in fractional calculus. In Section 3 Bernstein polynomials are introduced and 
their properties are shown in several lemmas. In Section 4, a FOCP with time 
delay will be solved using BPs. Section 5 contains some numerical examples. 
Finally Section 6 consists of a brief conclusion. 


2 Some preliminaries in fractional calculus 


Definition 2.1. A real function f(t), t > 0, is said to be in the space C,, 
u ER, if there exists a real number p > ys such that f(t) = t?fi(t), where 
fi(t) € C[0,+00) and it is said to be in the space Ci? iff f™ € C,, for 
meEN. 

Definition 2.2. The Riemann-Liouville fractional integral operator of order 
a > 0 of a function f € Cy, uw > 1, 1s defined as: 


olf f(t) = pty Jo @-— 7) “Far, 
olf f(t) = F(t). 


Definition 2.3. The fractional derivative of f(t) in the Caputo sense is 
defined as follows: 


(5) 


1 ; d” 
cpa ae a | mm 
6D f(t) a (t —T) gant) n-l<a<n, nEN, fec™ 


(6) 


In [15], the following properties for f € C,, and p > —1 have been proved 


1. of@t® = ret, keNu{o}, +t>0, 
2. §6Drolf f(t) = f@, 


aqlepDe iOS iO= Oe. SU, 


AC EDI FG) = ol EDA @, “ae 0. 
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3 Properties of Bernstein polynomials 


The Bernstein polynomial of degree n over the interval [a,b] is defined as 


follows: ; ; 
t-—a n\ (t—a\'(/b—-t\"" 
Bin (S=3) = CG) (Ga) (ia) 


so, within the interval [0, 1] we have 


Define ®,,(t) = [Bom(t) Bim(t) --: Bmm(t)|*. To consider the vector 
®,,(t — 7) ( 7 is the given delay ) in terms of ®,,,(t), we state the follow- 
ing lemmas. 


Lemma 3.1. We can write ®,,(t) = AT(t), where A = (Gao is 
an upper triangular (m+ 1) x (m+1) matrix with entry 


—1)I-t(™M (™) BG 
fag 45 ) es i,j =0,1,--- 


Proof. [4]. 


Lemma 3.2. For each given constant delay n > 0, ®n(t — 7) = 26,,(t), 
where Q is an (m+1) x (m+1) matrix in terms of 1. 


Proof. According to Lemma 3.1 we have 
®,,(¢— 7) = AT n»(t — 0). 
But, the right hand side of the above equation can be written as 


1 
t-— t 
2 2 

AT (t—n) =A| &-7) | =avl]* | =AvT, (0), 
(t—)™ ™ 


where 
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1 0 0 --- 0 

—n 1 0 --- 0 

a= n? —2n Li es830 
Speen ss A 


By Lemma 3.1, Tim (t) = A7+®(t), thus 


&(t — 7) = AVA~',,(t) = 06,,,(t). (7) 


Lemma 3.3. Let L?[0,1] be a Hilbert space with inner product (f,g) = 
i, f(t)g(t)dt and y € L?[0,1]. Then one can find the unique vector C = 


[Co C1 + °° Cm)" such that 
y(t) © Yo Bim (t) = 07 Bn() (8) 
i=0 


Proof. [12]. 


In Lemma 3.3 we have C? = Q-1(y, ®,,) such that 


1 
(y,®m) = Fi y(t) ®m(t)de = [(y, Bom) (y, Bim) «++ (Ys Bmm)]" 
and each entry of the matrix Q = (Q:41,341)4 50 is defined as follows: 


CG) 


(2m +1) (7) 


1 
Qisrger = | Bim ()Bim(t)ar = 
0 


Since the set {Bo m(t), Bim(t),-:: ,Bm,m(t)} forms a basis for the vector 
space of polynomials of real coefficients and degree no more than m [7, 16], 
a polynomial of degree m can be expanded in terms of a linear combination 
of Bim(t), @ = 0,1,---+ ,m) as follows 


moreover we have 


Lemma 3.4. Derivatives of P,(f) = pea, f(Z)Bjn(t) of any order converge 
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to corresponding derivatives of f. So if f € C*[0,1] then 
linea Palf))™ = f™, 


uniformly on [0, 1). 


Proof. [8]. 


4 Fractional optimal control problem with delays in 
control and state 


Consider fractional delay control system (2). For each 0 <i < r, one can 
apply the Riemann-Liouville fractional integral of to both sides of that 
equation 


aj(t) — 2; (0) = UF _ olf {ai,j (1); (t)} + Dperole {bin (tux (t) }+ 
Ui a0le { (aa)a,5 (t)a5(t — m)} + Se 102 {(ba)i,n (tur (t — bap 
9 


Assume that 2;(t) » X7®,,(t) (1 <i<r) and uz(t) UP ®»(t) (1 < k < s) 
where the entries X; = [X;(0)---Xj(m)]? and Ux = [Ux(0)---Ux(m)]? are 
respectively the coeffitients of x;(¢) and u,(t) in approximating them by Bern- 
stein polynomials of degree m just like (8). Moreover, the Bernstein approxi- 
mated coefficients vectors of functions a;,;(t), b:,4(t), (aa)i,j(t) and (ba):,4(t) 
can be achieved by using equation (8). We denote the approximated vector 
coefficients of these functions respectively by (A® 3)¢ (m+41)x 1 (BY) ntt)xds 
(Ay OY aise and (B PY eet hi 

By substituting the so called approximated vectors and matrices in (1), one 
can find the following equations: 


XP Om(t) — Bio = UVj-rolP{((A)7 Om (1) (AF Om(t))* 
+E parol (Bei) Om (t)) (Ug Om ( 
+Eja10Lf {((AG?)7 Om (O))(XF Om(t — my) 
+E olf {(B Bi a )7 Om ())(UZ®, ( 


Moreover, from Lemma 3.2 there exist (m+ 1) x (m+ 1) matrices Q),Q2 
where ®,,(¢ — 71) = 1 ®»(t) and ©®,,(t — 42) = Q2Pm(t), while 


0, = AVAT, 
= AW’ A7}, 


and U,V’ are obtained respectively in terms of 7 and 1. 
As it was shown in [4], foreach 1 < i,j <rand1<k< s, the (m+1)x(m+1) 
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matrices A‘ : Bik, Ais and Bee can be calculated such that: 

(A)7S,, ()®F, (t) = OF, QA, 

(BEF)? ®,, (HOF (t) = OF (OB, 


(Ai7)P Om (OF (t) = OF (t) AU, 


(Bi)? ®,,(t) 82 (t) = OT (t) BY". 


Therefore, by replacing the above equalities, (2) can be rewritten as follows: 


XFOn(t) — 4,0 = Dhna (ol (Q)(AMIX;) + Dh, (LPB, ())(BUU a) + 
Dyna (olf f(t) (AW OT X5) + Uh (olf &, (1))(By? OF Up), 


or 


XP @m(t) — 24,0 = Eh (AbIX;)T (olf Bm (t)) + Bha (Be*U a) (oP Bia(t))+ 
Dh (AWOL X;)7 (of Om (t)) + Dh (By OFUg)? (olf in| 


where i= 1,--- ,r. 


One can approximate oI?'®,,(t) by Ig, x ®y,(t), where I, is an (m+1) x 
(m+ 1) matrix called the operational matrix of Riemann-Liouville fractional 
integral. 

Infact, from Lemma 3.1, ®,,(t) = AT,(t), so 


ol? P(t) = A ol? Tm(t) = A [o[21 olSt --- ol Pt™]’, 


where o/?t? = rot. Therefore, 


of T(t) = BF, (12) 


where © = (Si41j41) and T = (Tj41) are respectively (m+ 1) x (m+1) and 
(m+ 1) x 1 matrices, which are defined as follows: 


2 D(j+1) = j 
Mita g+1 = ae , 1,9 =0,--+,m 


and 


Qs Sr 2S 


rss 5m. 


Also, from Lemma 3.3, since t’t? € L?((0,1]) for each integer i (0 < i < m), 
one can find the (m+ 1) x 1 vector P; such that 


ion? co, PPO. (4), (13) 
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where P; = Q-'{t't*,®,,(t)) and the entries of P, = (t't*,®m(t)) = 


[P;.0 Pia +++ Pim]? can be attained as 


é te gg Tatjtatl 
Pax f t+? BR. (t)dt = (i+j+atl) Pte an 
0 


— TG +m+at 2)’ 


Now if P is an (m+1) x (m+1) matrix of the form [Po Pi --- Pm], then from 
(12) and (13) we have 


ol, (t) ¥ AXPTS,,,(t), (14) 


therefore, Ig = A™PT is the aforementioned operational matrix of Riemann- 
Liouville fractional integral oI}. 


Hence, by replacing o[?®,,,(t) from (14) into (4) and writing x; in terms 
of BPs of degree m, equation (4) can be written as the following 


XP ®,,(t) — XZpEn(t) = Dt, (AMI XG)? Tash) De_,(BYRU,)? I. ®m(t) 
+054 (AG? OP X;)7 Ig Py (t) + He (BOL, IP, (), 
(15) 
where 
XT, = [Xi,0(0),--- ,Xi,o(m)]? 


is the known Bernstein approximated coefficents vector of x;,9 that can be 
computed using (8). By equalling the coefficents of ®,,,(¢) from both sides of 
(5), we found that 


AP exis Deg X? AG) 1, +E BM) Ie 
+25 XF01(Aq’)* Ta + wpa Ug N2(By’)” Ta, 


(16) 
for i =1,--- ,r. Equations (8) can be written in compact form as follows: 
XP? =11+uU"'T, (17) 


where II and I are respectively 1 x (m+ 1) and (m+ 1) x (m+ 1) matrices 
that can be obtained by the following 


RaeX GAL ADR: 


and 2 a : 
T=(B+ Baldi — (A+ Aaa), 


and Im+1 is the (m+ 1) x (m+ 1) identity matrix. 
Moreover, by applying the approximations x(t) © (XT)1.-¢m+1)®m(t) 


and u(t) * (U*)ixs(m+1)@m(t) where X7 = [X7,---,X7] and UT = 
[U7 ,--- ,UZ], the cost functional (1) can be approximated as bellow 
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J=i Jy {27 )QOzx(t) +ul (t)R(t)u(t)}dt 
= HDX Hnl VQ? Bul OBROAY +O Onl N(R Bn (OPH) ae 
18 
where Q = [Q;,;] and R = [R;,;] that Q;,;, Rij; are the (m+ 1) x 1 vec- 
tors of Bernstein coefficents in approximating Q,,;(t) and R;,;(t) respectively. 
Therefore, 


Tah [xT a (HBL OXY + UP Hm(N OL ORV)” Ht 
or 
oe 5 [ (XT mNXTOTHn(O) + (UTOm(t))(UT RT On (t)) }dt, (19) 
0 


where Q = [Q;,;] and R = [R;,]. Also Q;,; and Rj; are (m+1) x (m+1) 
matrices that can be calculated from 


(Q*4)7 Om (1) BF, (t) = BF, (NQ*9, 
(RY)? On (OR) = OR (QRY. 


Let Zi; = H® Qi; and Wi; = H@R,;, where ® is the Kronecker 
product and H = [Ajj)¢m+41)x(m+1) and each entry H;,; is defined by 


1 
Hi,5 = | Big Benlt dt, 
0 
then (19) can be rewritten in compact form as: 


Jz 5 {(X7ZX) + (UTWU)}, (20) 


where Z = [z;,;] and W = [w,,;]. 
From (17) we know that X7 = II + U7T, so the necessary condition that U 
minimizes (20) and satisfy (17) is that 


OJ 

— = xXTZIT+UTW =0 

aU ‘is 
so 


ae Gale A a a (21) 


The above equation gives the optimal feedback control and by replacing (21) 
in (17), we can easily find the optimal state as well. 


We need to mention that, since the Bernstein coefficients of positive func- 
tions in L?[0, 1] are positive [7] and it was assumed that R(t) is positive def- 
inite, then R;,; is a positive vector. Also because B;m(t) > 0 for t € (0,1), 
it’s clear that 
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OF (t)Rij = RE; Om(HOR(t)>0, te [0,1], 


therefore Ri; and as the result W;,; = H @ R;,; are positive definite and 
consequently invertible matrices. 


5 Convergence of the method 


In this section, we show the convergence of the presented method discussed 
in this article. First we prove the following lemma. 


Lemma 5.1. Let X7,,(t) = 0 Xj Bjm(t) be the Bernstein poly- 
nomial of order m that approximates the function x(t) € L?[0,1]. Then 
ol*(X7S,,,(t)), tends to oI*x(t) as m tends to infinity. 


Proof. By Lemma 3.3 we have 


lisse DGB a). (22) 
j=0 


Since B;(t) is a continuous function, we have 


* i-0 Xj Bi m(7) mn * Bym(r) 
limm of Jee dt = limm-soo x; f dr. 
a [Fatal ay tin Sf Path 


By (22) and from Definition 2.2, we obtain 


t m 

| ty ~dt =T (a) limmsoo S° X; olf Bim(t), 
a a 

or 


olfx(t) = limmsoo )) Xj olf Bjm(t) = littm+ooX oh“Bm(t). (23) 
j=0 


In (14), I, = ASP? where the i-th column of P? is the Bernstein approx- 
imated coefficients of t‘t® for i = 0,--- ,m. Now, regarding the convergence 
of the Bernstein approximation of every functions in L?([0,1]), one can write 
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A 
» 5-0 P,0Bj,n(t) site 
= PjiB n(t) ~ +a 
littnsc0P 7 @n(t) =litnsco | Spa er" ||, 
ve i=0 PjmBjn(t) pmta 


therefore, limn—+ooADP? ©, (t) = A> litttinseeP? O,(t) = AST: or as ex- 
plained in (12) and (13) 


limn—+cola Pn (t) = ols Pm (t). (24) 
From (23) and (24) we reach 


olea(t) = litimsoX? limn+colaPn(t). 


Given n > m will complete the proof. 


Theorem 5.1. The approximated solutions z(t) = X7®,,(t) and u(t) = 
UT®,,(t) in which (X,U) is achieved from (17) and (21), converge to the 
optimal solutions «*(t) and u*(t) as the degree of the Bernstein polynomials 
tend to infinity. 


Proof. Suppose W,, is the set of all (U7,X7)®,,(-) where X,U € R™* 
and satisfy (17), also W is the set of all (u(-),2(-)) satisfy (2) and (38). 
Let U be the optimal solution of (20) where obtained from (21) and X be 
the solution of (17) obtained by replacing U in eqation (17). Therefore 
(UT, X7)®,,(-) € Wm. By the convergence property of Bernstein polyno- 
mials, for (U7, X7)®,,(-), there exists a unique pair of functions (a(-), Z(-)) 
such that 
(OT, X7)®,,(-) — (a(-),(-)) as m— oo. 


Now according to Lemma 5.1 it is clear that (u(-),Z(-)) € W. Moreover as 
m —> oo, then J(U'S,,, X7®,,) —> J where J is the value of cost function 
(1) corresponding to the feasible solution (t(-), Z(-)). Now, since 


WiC: CWm C Wing CCW, 


consequently 


Infw,J1 SS InfwIm oa Ln fWing1Im41 St OS Infw J. 


Let J*, = Infw,,Jm, so J*, = J(UT®,,, X7©,,). Furthermore, the sequence 
{J} is nonincreasing and bounded bellow which converges to a number 
J > InfwJ. We want to show that J = limm+5.J*, = InfwJ. Given 
é > 0, let (u(-),2(-)) be an element in W such that 


J(u,2) <InfwJ +e, (25) 
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by the definition of infimum, such (u(-), (-)) € W exists. 
Since J(u,2) is continuous, for this value of ¢, there exists N(¢) so that if 
m> N(e), 

|J(u, x) — J(UT®m, X7Om)| <eE, (26) 


Now if m > N(e), then using (25) and (26) gives 
JUF Oy; X7 On) < J(u,2) +e < InfwJ + 2¢, 
on the other hand 
InfwJ < J*, = Infw.,Jm < J(UT Em, X7 Bm), 


sO 
Infwd < Jp, < InfwJ + 2¢, 


or 
O< I Hind Soe 


where ¢ is chosen arbitrary. Thus 


J =limmsoJ7, = Infw J. 


6 Numerical examples 


In this section we give some numerical examples and apply the method 
presented in Section 4 for solving them. Our examples are solved using 
Matlab2011a on an Intel Core i5-430M processor with 4 GB of DDR3 Mem- 
ory. These test problems demonstrate the validity and efficiency of this tech- 
nique. 


Example 6.1. Consider the following delay fractional optimal control prob- 
lem in which 0 <a <1, 


min J =} f) [x(t) + du2(a)]de, 

st §D@a(t) = —a(t) + a(t- 3) + ud) — fu(t- 3), O<t<1, 
a(t) = 1; -i<t<0, 
u(t) =0, -l<t<o 


For a = 1, this problem has been numerically solved by applying hybrid 
functions based on Legendre polynomials in [19] and the objective value 
I = 0.3731 has been achieved. Whilst, in the presented method the solu- 
tion has the objective value J* = 0.3956 for a = 1 and m = 6. Thus, our 
results with m = 6 are in good agreement with the results demonstrated in 
[19] for a = 1. In addition, by varying the value of a we can obtain the op- 
timal control u(-) and trajectory function x(-) which are shown respectively 
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Figure 1: Approximate solution of u(.) for a = 1, 0.999, 0.99 in Example 6.1 


Figure 2: Approximate solution of a(.) for a = 1, 0.999, 0.99 in Example 6.1 
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Table 1: The objective value and the end point of trajectory for a = 1,0.999,0.99 in 
Example 6.1 


a | objective value _ end-point 

1 0.3956 0.6775 
0.999 0.3283 0.6443 
0.99 0.2907 0.6249 


Table 2: The objective value and the end points of trajectories for a = 1,0.9,0.8 in 
Example 6.2 


a | objective value end points 

1 0.7245 —0.4691 , — 0.0113 
0.9 1.0291 —0.6477 , 0.3202 
0.8 0.7299 —0.4324 , 0.4674 


for some values of a in Fig.1 and Fig.2. Moreover, for these values of a the 
objective values and the end points of optimal trajectory are shown in Table 
1. 


Example 6.2. Consider the following two-dimensional DFOCP in which 
0<a<l, 


min J = 4 f{[xr(t) x9(t)| E 2| [a y(t) xo(t)|” + (t2 + 1)u2(t)} dt, 


ote (8)- [Fo] ERREH] ef Jon [Jebs oes 


[vi(t) vo(t)] = [1 Y, — 


u(t) =1, -— 


dle 
IA 
IA 
=) 


<t<0. 


Ale 


This problem for w = 1 has been studied in [19], where the obtained approx- 
imated cost function is J = 1.5622. Using the presented method for a = 1 
and m = 6, gives the approximated cost function as J* = 0.7245. So we 
achieved satisfactory numerical results in comparison with what have been 
obtained in [19] for a = 1. Also by varying the value of a the obtained 
control and trajectories functions are shown respectively in Fig.3, Fig.4 and 
Fig.5. Moreover, for these values of a the objective values and the end points 
of optimal trajectories are shown in Table 2. 
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Figure 3: Approximate solution of u(.) for a = 1,0.9,0.8 in Example 6.2 


Figure 4: Approximate solution of 21(.) for a = 1,0.9,0.8 in Example 6.2 
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Figure 5: Approximate solution of x2(.) for a = 1,0.9,0.8 in Example 6.2 


7 Conclusion 


In this paper, we peresent a new method of using Bernstein polynomials 
for solving DFOCP’s. We approximate the objective function and find a 
feed back control which minimizes the cost function. Then by replacing the 
optimal control in the constraints, we get an algabric system which can be 
solved in terms of the approximate coefficents of trajectory. The convergence 
of the method is extensively discussed and some test problems are included 
to show the efficiency of this very easy to use and accurate method. 
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